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Description 

[0001} This invention relates in general to medical imaging, and more particular to an improved segmentation 
method for 3-dimensional ultrasound (3-D US). 

5 [0002] The utility of ultrasonography in the diagnosis and assessment of carotid disease is well established. 
Because of its non-invasive nature and continuing improvements in image quality and Doppler information, ultrasonog- 
raphy is becoming increasingly popular in such applications. In 1991, the NASCET results demonstrated that an angi- 
ographic stenosis of 270% selected a group of patients that benefited from carotid endarterectomy, but no equivalent 
Doppler measurement could satisfactorily select this group {see R.N. Rankin, A.J. Fox, K. Thorpe, and NASCET col- 

10 laborators, "Carotid ultrasound: correlation with angiography in a multicenter trial," Presented to the American Society 
of Neuroradiology, St. Louis: 31 May - 5 June (1 992)). Although the role of ultrasound as a screening tool is well estab- 
lished, its role as the only definitive diagnostic test in assessing risk for stroke before surgery is very controversial and 
subject to heated debate. Nevertheless, there is agreement that factors related to the flexibility of ultrasonography and 
its high machine and operator dependence have contributed to disappointing results in some tests and trials. Nonstand- 

15 ardized techniques and improper choice of equipment yield variable results, particularly in multi-centre trials. In addi- 
tion, measurement choices have in some circumstances been other than optimal. Flow-velocity-based measurements 
of stenosis severity are subject to a great deal of variability due to measurement location and angle uncertainty, as well 
as variations in operator skill. These variabilities and inaccuracies are further increased because the measurement of 
a single velocity vector component at one or two locations in a vessel constitutes only an indirect measure of the risk of 

20 stroke. 

[0003] It has been speculated that 3-D US with real-time visualization of plaque and its surface, as well as 3-D 
measurements of the stenosis and the actual atheroma volume would reduce the variability of carotid disease assess- 
ment and improve diagnosis (see (1) T.S. Hatsukami, B.D. Thackray, and J.F. Primozich, "Echolucent regions in carotid 
plaque: Preliminary analysis comparing three-dimensional histologic reconstructions to sonographic findings," Ultra- 

25 sound Med & Biol 20, 743-749 (1994); (2) PH Arbeille. C. Desombre, B. Aesh. M. Philippot, and F. Lapierre, "Quantifi- 
cation and assessment of carotid artery lesions: degree of stenosis and plaque volume," J Clin Ultrasound 123, 1 13- 
1 24 (1995); (3) W. Steinke. and M. Hennerici. Three-dimensional ultrasound imaging of carotid artery plaques," Journal 
of Cardiovascular Technology 8, 15-22 (1989); and (4) D.D. McPherson, "Three-dimensional arterial imaging," Scientific 
American Science & Medicine 22-31. (March/April 1996)). 

30 [0004] In addition, 3-D imaging of plaque may also allow quantitative monitoring of plaque development (such as 
changes in volume and morphology), provide important information about the natural history of atheroma growth, and 
help in the identification of those plaques which represent a risk of giving rise to stroke. 

[0005] With a 3-D ultrasound image of the carotid arteries, important information about the vessels, which previ- 
ously was difficult or impossible to obtain accurately and reliably, can now be ascertained after examination of a patient 

35 In the measurement of the stenosis, the diagnostician can step through the 3-D image, one slice at a time, and outline 
the edges of the vessel wall. Counting the number of pixels enclosed within the trace and multiplying by the area of a 
pixel gives the cross-sectional area of the vessel. However, this manual process is highly labour intensive and operator 
dependent. It has been recognized that a preferable approach is to utilize automatic segmenting of the vessel. 
[0006] Much work has been reported on the semi-automated segmentation of ultrasound images. This includes 

40 work on left ventricle boundary detection in echocardiographic images, ovarian follicle extraction, and intravascular 
ultrasound segmentation. While many of these studies have produced useful results, nearly all require the user to spec- 
ify an initial contour, giving rise to inconsistent and mixed results from user to user. Additionally, most are computation- 
ally intensive and cannot be practically applied to a large 3-D image. 
[0007] The invention is set out in the claims. 

45 [0008] According to the present invention, an automated segmentation method is provided for three-dimensional 
vascular ultrasound images. The method includes two steps: an automated initial contour identification, followed by 
application of a geometrically deformable model (GDM). The formation of the initial contours involves the input of a sin- 
gle seed point by the user, and has been found to be insensitive to the placement of the seed within a structure. The 
GDM minimizes contour energy, providing a smoothed final result, with only three simple parameters being required as 

so easily selectable input values. The method according to the present invention is fast (capable of performing segmenta- 
tion on a 336x352x200 volume in 25 seconds when running on a 100 MHz 9500 Power Macintosh prototype) and 
involves minimal user interaction and minimal processing. The inventive method addresses the prior art problem of var- 
iability created trough user interaction, particularly due to the definition of initial contours and the choice of threshold 
parameters. 

ss [0009] Embodiments of the invention will now be described by way of example, with reference to the drawings of 
which: 

Figure 1 is a schematic representation of the Geometrically Deformed Model (GDM) of a closed contour formed 
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from a series of vertices connected by vectors. 

Figure 2 is a schematic representation of a hypothetical vessel boundary, with search rays extending radially from 
a single seed point for the purpose of edge selection during an automated initial contour identification step accord- 
ing to the method of the present invention. 
5 Figure 3 is a graph showing an intensity profile and the corresponding Threshold Tree with five threshold levels for 
the hypothetical vessel boundary of Figure 2. 

Figure 4 is a schematic GDM curvature representation (a) at a vertex defined as the subtraction of two adjacent 
vectors, and (b) at three angles, for minimizing contour energy following the initial contour identification step 
according to the present invention. 
10 Figure 5 is a schematic representation showing the tangent and radial vectors at a vertex, for calculating external 
force used in the GDM. 

Figure 6 shows an example of the results of applying the automated segmentation method according to the present 
invention, wherein Figure 6(a) is an image of a single vessel slice showing the placement of 1 5 seed points, Figure 
6(b) is a plot of the initial contour, and Figure 6(c) is a plot of the final contour. 
15 Figure 7 shows the effect of the GDM application step on contour energy in a 60% stenosed phantom vessel at a 
point above the bifurcation, wherein Figure 7(a) is a graph showing the total energy, and the internal and external 
components as a function of iteration number, Figure 7(b) shows the initial contour prior to application of the GDM 
step, and Figure 7(c) shows the final contour after the last iteration. 

Figure 8 shows the segmentation results according to the method of the present invention for a 60% stenosis phan- 
20 torn at the common carotid (Figure 8(a)), at the bifurcation (Figure 8(b)), and above the bifurcation at the internal 
and external carotids (Figure 8(c)). 

Figure 9 is a multi-planar reformatted 3-0 image of the 60% stenosed phantom of Figures 7 and 8 showing three 
cut planes. 

Figure 10 shows lumen cross-sectional area as a function of distance along each branch for both the numerical 
25 models used in casting the phantoms used in the examples of Figures 7 - 9 and the segmented results for (a) 30% 
stenosis; (b) 60%; and (c) 70%. 

Figure 1 1 shows the effects of GDM resolution on the contours at a point above the bifurcation in the 60% stenosed 
phantom, wherein Figure 1 1(a) shows an initial contour before application of the GDM, Figure 1 1(b) shows the final 
contour for a vertex spacing of r = 8 pixels; Figure 1 1(c) shows the final contour for a vertex spacing of r = 1 5 pixels; 
30 and Figure 1 1(d) shows the final contour for a vertex spacing of r = 25 pixels. 

Figure 12 shows the results of the segmentation method according to tha present invention on patient carotids 
acquired with (a) a mechanical scanner and (b) a freehand system using a magnetic positioning device. 

[0010] Before describing the preferred embodiment and best mode of the invention, prior art methodologies of 
35 active contour models and GDM are presented briefly herein since aspects of these methodologies are incorporated 
into the automated segmentation method according to the present invention. 

[001 1] Segmentation involves separating an image into its constituent parts. Many techniques are known in the art, 
some of which have been described briefly above. However, nearly all such prior art techniques utilize a limited number 
of basic methods including gradient operators, intensity thresholding, template matching, and region-based analysis 
40 such as region growing (see RC Gonzalez and P. Wintz, Digital Image Processing, (Addison-Wesley, Reading, MA, 
1987)). These techniques rely either on changes in intensity across a boundary (local operators), or on the similarities 
of the region within a boundary (regional operators). In medical imaging, a wide diversity of structures may be encoun- 
tered with little regularity, even within the same class of structures. Although organ images of different individuals are 
similar, individual and pathological differences increase the diversity. This makes it difficult to find a general represen- 
ts tation for the structures using fixed features or templates. Traditional edge detectors are even less useful in ultrasound 
imaging where speckle, shadowing, and other intrinsic noise and image artifacts degrade the image quality. Detectors 
that rely on local operators such as gradient, are highly sensitive to noise and artifacts. Region-based operators are less 
sensitive. However, they can be biased by discontinuities in the data or by gradual changes in contrast, such as from 
shadowing. Each approach relies on different properties of the image, but no technique alone is satisfactory in bound- 
so ary identification and segmentation. When the properties of local and regional detectors are combined, such as in the 
active contour and the geometrically deformable models, the quality of the analysis can be improved. 

A. Active Contour Models 

55 [0012] The "snake" active contour model was first introduced by Kass et al (see M Kass, A. Witkin, and D. Terzo- 
poulos, "Snakes: active contour model," International Journal of Computer Vision 1, 4:321-33 1 (1988)) thereafter 
refined by others. The "snake" active contour model is one of the most widely used segmentation methods presently in 
existence. This method is well suited for applications where the object of interest has a natural, continuous, and smooth 
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contour, as often found in medical images. Generally speaking, the models are semi-automated, requiring the user to 
provide an initial close estimate of the contour within the true contour boundaries. This initial contour can be interpreted 
to contain the regional information. The curve evolves through a series of deformations as influenced by the attraction 
of external forces of nearby image features such as edges and lines, while simultaneously maintaining smoothness 
5 constraints on the evolving shape of the contour, which is controlled by the internal forces within the contour. The inter- 
nal forces are a function only of the contour's curvature, and the goal is to minimize the overall curvature while still main- 
taining a degree of rigidity. 

[0013] The first generation of snake models (and the most common) uses the normalized arc length s to parame- 
terize an initial curve as v(s)=(x(s),y(s)) , where s varies from 0 to 1 along the curve and x(s) and y(s) are the corre- 
10 sponding coordinate points along its length. The total force. F, oto) , acting on the contour is the sum of the internal and 
external forces: 

1 

F total = f F M( v ( s ))* F e X ,( v ( s )))< is 



[0014] The internal force, F lnt is usually defined by a combination of the first and second derivatives along the 
curve. The first derivative ensures that the contour remains continuous and flexible (like a rubber band), while the sec- 

20 ond derivative imposes a degree of stiffness, keeping the contour smooth and preventing it from collapsing upon itself 
The external force, F 6xt is a function of the image properties, usually defined as the local gradient at each point along 
the curve. The curve's evolution is guided by the solution that minimizes the total acting forces. The solution is obtained 
iteratively, and is generally solved through variational calculus methods. The contour slowly approaches the true bound- 
ary, with each point along the contour adjusted during each iteration. The optimal solution is mathematically intensive 

25 and computationally slow. 

[001 5] Other drawbacks to prior art active contour models include their sensitivity to the initial contour, as well as to 
the numerous parameters, which must be chosen by the user. If the initial contour is too far from the object edges, the 
contour will not be attracted by the distant gradients. 

30 B. Geometrically Deformable Models 

[0016] Geometrically Deformable Models (GOMs), introduced by Miller et al.. are a different class of active contour 
models which avoid some of these complexities (see J.V. Miller, D.E. Breen, W.E. Lorensen, R.M. O'Bara, and Wozny 
MJ, "Geometrically deformed models: a method for enacting closed geometric models from volume data," Computer 

35 Graphics 25, 4:217-226 (1991)). 

[001 7] While GDMs share many similarities with the "snake" models discussed above, they differ fundamentally in 
how they define the contour and, accordingly, how the deformation occurs. The "snake" models are basically a series 
of connected spline segments. During deformation, optimization is performed along the entire spline and all points are 
adjusted. In GDMs, the contour is defined as a series of vertices (Figure 1) which are in turn connected by edges. This 

ao is in essence a discrete model, with the resolution determined by the length of the edges between vertices. The more 
simplified discrete structure of GDMs allows the internal and external forces to be calculated independently at each ver- 
tex using local information. The deformations are applied at the vertices only, with the goal of minimizing the forces act- 
ing upon each vertex. 

[0018] The GDM's force function, in particular the internal force component also differs from that of the snake 
as model. For the GDM, the curvature at a vertex is approximated by the magnitude of the vector difference between the 
edges joining at that vertex. The simplicity of the underlying structure of GDMs allows for an easier and much faster 
convergence to a final contour than prior art "snake" models. This simplicity also translates to a reduction in the number 
of user-specified parameters to just three: the internal and external force weightings, and the between-vertex spacings. 
[0019] According to the method of the present invention, aspects of both the local and regional detectors found in 
50 the simpler and more efficient GDM technique have been adopted. In order to address the large source of variability in 
the segmentation arising in the prior art from the requirement for a user-defined initial contour, the method of the 
present invention uses an automated contour initialization and require the specification of only a single seed point in the 
approximate center of the vessel lumen in any slice of the vessel. This automated initial contour selection is particularly 
well suited to the subsequent GDM application step which is sensitive to the placement of this initial contour (i.e. the 
55 external forces, which pull the contour to the true boundaries, are local in nature and can exert their forces only within 
a limited region). 

[0020] In the following description, the various stages of the inventive segmentation algorithm are presented, 
including initialization and deformation, followed by a discussion of experimental results. 
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1. Initialization 

[0021] For the case of transverse images of the carotid arteries or any other fluid-filled structures, the contours will 
be closed or nearly closed. As such, a two-step method is provided according to this invention for automatic initial con- 
5 tour generation. 

a. Edge Selection 

[0022] A collection of straight line rays 5° apart are extended radially outward from the seed as shown in Figure 2. 
w A preprocessing operation is performed along the ray using a 3x3x3 kernel Specifically, at each voxel along the ray, 
l(i,j,k), a local median, mean, and variance operation is performed on the 3x3x3 kernel, generating the values md(ij.k), 
mn(i,j,k), and v(i,j,k), respectively. If the variance between the voxel's intensity and the local mean is 50% greater than 
the avenge local variance, it is replaced by the median value of the kernel according to the following two equations: 

var(/,/,fc) 

ifr(i,j,k)>o.5 

20 then l{i,j,k) = md(i,j,k) 

[0023] The intensity profile of the resultant ray is then searched for the boundary of the vessel using threshold anal- 
ysis. The profile is evaluated at multiple thresholds, resulting in a multilevel binary Threshold tree" as shown schemat- 
ically in Figure 3 for the typical profile of an edge. A "0" indicates a voxel intensity below the threshold value and a "1" 
25 above. 

[0024] To avoid variability due to user-defined parameters, thresholds are automatically selected for each ray, mak- 
ing them adaptive to the image properties. A simple search is conducted for the maximum intensity along the ray, with 
the condition it any maximum conforms to a trend and is not an isolated noise occurrence, for calculating the five thresh- 
olds of the tree in Figure 3. In calculating these five thresholds, it is assumed that the voxel with maximum intensity, l m 
30 is likely a part of the boundary edge along this ray, and therefore the thresholds are chosen so as to have this voxel 
occupy at least four levels of the threshold tree. Using the intensity l 0 at the origin of the ray, along with the previously 
found l m , the thresholds T, are specified as follows: 

Al = (l m -l 0 )/4 

35 

T,*l o +<rA0i = 1,2,3.4,5 

[0025] The information in the "threshold tree" is used to create a list of potential candidates for the vessel wall 
boundary. If an analysis of the resulting candidate edge points reveals that there are too many candidates, then Al is 
40 increased using Al = (I m - 1 D )/3 , and the threshold tree is reevaluated. 

[0026] From the information in the tree, three descriptors of a candidate edge are determined. These are: start 
position, length, and number of levels occupied. For the tree in Figure 3, the result is as shown in Table I. 



Table I 



Results of profile analysis for edge in figure 3 


Candidate 


Position 


Number of Levels 


Edge Length 


1 


4 


1 


1 


2 


6 


2 


1 


3 


9 


5 


7+ 



55 [0027] Using these descriptors, a set of four criteria are derived for rating the candidates. These are: edge position 
relative to that of possible candidates from other adjacent rays within the slice, the edge strength (or the number of lev- 
els occupied), the edge length, and the variance in the voxels along the ray between the origin and the candidate edge. 
Using these criteria, an initial boundary point is selected as described, and the points are then further refined as 
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described below under the discussion of Contour Formation. 

[0028] More particularly, within a set of candidates in a ray. each criterion is normalized across ail candidates to 
give a value between 0 and 1 . This ensures that at least one candidate receives the maximum value in each of the cri- 



> nm(£;.,.,) 



where i=1 to 4, representing the four criteria indicated above. 
[0029] For an edge, the sum total would be given by: 



where Wj is the relative weighting of criterion i. Initially, w, = 1 for all criteria. If the point total for one of the edges, Sj, is 
20 significantly greater than any of the others, it is automatically selected as the initial boundary point: 

if Sj > 1 + (maxS,,) , k = 1 ton, k*j 
then j = initial edge 

25 [0030] If no dominant edge is present, the relative weightings of the criteria are readjusted. Difficulty in selecting an 
edge occurs most often in areas of strong shadowing and noise artifacts. In these areas it is more useful to place the 
greatest emphasis on proximity (criterion 1), followed by length (2), strength (3), and variance (4). Therefore, initial 
weight adjustments are chosen to reflect this priority, and are set as follows: 

30 w, = 1.0 



[0031] Subsequent weight adjustments are set as w ( = Wj - 0.1 , i = 2 to 4, with the proximity criterion, w h remaining 
35 at 1.0, throughout. In this way, an attempt is made to include as many of the edge criteria as possible in selecting the 
boundary, while slowly shifting the main emphasis towards the proximity criterion and the continuity of the contour. 

b. Contour Formation 

40 [0032] The proximity, or continuity, criterion constrains adjacent rays to select boundary points that are continuous 
within a slice. To extend this continuity constraint to three dimensions, a second pass is made through the data and the 
boundary placement is finalized. Each previously selected boundary location is compared with that of the candidates in . 
the two adjacent rays in the current plane and the three in each of the adjacent planes, giving eight neighboring bound- 
ary points. The selected point is then compared to its eight neighbors, and if it is outside the volume defined by them., 

45 then other candidates are examined. Otherwise, the point is considered to be the best boundary location and is saved 
as a point along the forming contour. 

[0033] At the end of this process, a set of automatically generated contours is provided consisting of discrete points 
connected by edges as shown in Figure 1. This is precisely the configuration needed by the GDMs. 

so 2. Deformation 

[0034] According to the segmentation method of the present invention, the well known deformation technique 
described by Lobregt et al.. has been adopted (see S. Lobregt and M.A. Viergever, "A discrete dynamic contour model," 
IEEE Transactions on Medical Imaging 14, 12-24 (March 1995)). The vertices illustrated in Figure 1 represent the i m 
55 vertex with cartesian coordinates (v ix , Vjy). The vectors dj connect adjacent vertices, and are labeled in a counter-clock- 
wise direction. The length of the vectors d,, (||di||), connecting the vertices determines the local contour resolution. If 
the vector lengths are allowed to be too large, small changes in the image cannot be followed, but if they are too small, 
the contour may track small noise variations. The length of d| can change at each iteration during deformation, but is 
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kept within bounds of ± 1/3 of a user-specified resolution, r. by removing or inserting vertices as needed. If ||d,|| exceeds 
the maximum allowed resolution, a new vertex is inserted midway between the two vertices, and if ||d||| is less than the 
minimum allowed resolution, then one is removed and the other is adjusted to lie midway between its original position 
and that of the removed vertex. 

5 

a. Internal Force 

[00351 With a vector model, the internal force can be represented by the vertices' curvature. The curvature can be 
described very simply as the vector difference between the two edges sharing that vertex. Using the normalized vectors 
10 d,, and d w , Figure 4a illustrates how vector subtraction gives a curvature vector C|, that provides a measure of the angle 
between the two edges. As illustrated in Figure 4b, the sharper the local curvature, the longer this vector will be. It 
should be noted that the magnitude of the curvature vector, varies between zero and two. In order to minimize the 
internal contour force, the vertex v, must be moved in the direction of the curvature vector, making 

' 5 F* = c r 



[0036] In minimizing the curvature of the contour as a whole, areas of constant or slowly varying curvature should 
20 not be affected. To do otherwise would cause circular and polygonal shapes to continue to pull the vertices inward until 
they converge into a single point This problem is avoided by applying a symmetrical filter centered at each vertex, and 
making the internal force a function of the filtered curvature vectors, as set forth in the Lobregt et at paper. According to 
the preferred embodiment, the following filter parameters were used: 



Therefore, the internal force at vertex / would be 

30 

F, n( =(-VS||c M |H|c I .||-i/J||d w ||) jj ^ jj 



3S b. External Force 

[0037] To find the external force at each vertex, the gradient g, is first calculated using the Sobel operator and the 
result is then transformed according to an exponentially decaying function to give a scaled gradient G, 

[0038] This equation was chosen based on observations over a wide range of images, where it was noted a gradi- 
ent magnitude less than 15 rarely occurred in regions on the border or outside of the structure of interest. A factor of 

45 two was used to give the scaled gradient a magnitude range of zero to two, as with the curvature vector. 

[0039] The vector Gj at a vertex V| can be separated into two perpendicular component vectors, one tangential to 
the contour and the other in the radial direction. It should be noted that moving the vertex along the contour (i.e. in the 
tangential direction) does not change the contour shape. Therefore, only a force in the radial direction is required to be 
defined in order to move the vertex radially. In Figure 5, it will be noted that the tangent vector at vertex v, is represented 

so by the sum of the normalized distance vectors 

t.-d. + d,.,. 



55 If the tangent is rotated by - jt/2, it will point in the local radial direction: 
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[0040] Therefore, the curvature and the local radial direction vectors are equivalent, causing the internal and exter- 
10 nal forces to act along the same radial line, either in the same or in opposite directions. To calculate the external force, 
the transformed gradient is normalized along the radial direction, thus: 

c. Total Force and Deformation 

[0041] The total force acting upon a vertex is a weighted sum of the internal and external forces, thus: 

20 

P -B i"w 0 F tl +w h F l . 



where the two parameters w„ and w^, determine the relative weightings. If the force at a vertex is sufficiently small such 
that it experiences no translations for three consecutive iterations, that vertex is deemed "stable" and is locked into its 
position. According to the best mode as of the time of filing this application, the deformation process is stopped when 
at least 90% of the vertices have been locked into place, or a maximum 100 iterations of the algorithm have been car- 
30 ried out. This constraint is imposed in order to ensure a controlled exit for a case where the initialization step has pro- 
duced a poor contour, with little chance of the vertices settling at a stable boundary. 

3. Experimental Results 

35 [0042] The 3-D images used in testing the algorithm of the present invention were all obtained using a Sims 3-D 
Ultrasound imaging system (Life Imaging Systems, Inc., London, Canada) with a liner scanning mechanism, as et forth 
in commonly assigned issued U.S. Patent Nos. 5.562,095 and 5,454,371 , the contents of which are incorporated herein 
by reference. This mechanism consists of a mounting device to hold a transducer in place, and a motor system to trans- 
late the transducer across the surface in a linear manner. The ultrasound machine used in the experiments was an ATL 

40 Ultramark-9 (Advanced Technology Laboratories. Bothell. Washington) with an L 10-5 38mm linear transducer array. 
The mounting device was manually adjusted to be perpendicular to the scanning direction, thus making the ultrasound 
beam transverse and perpendicular to the vessel. The scanning speed was set at 3. 75 mm/sec and the acquisition rate 
was set at 15 frames/sec, giving a between-frame separation of 0.25mm. A total of 200 2-D images covering a length 
of 50mm were acquired in about 13 seconds. The 2-D images were reconstructed to form a 3-D volume which was 

45 viewed and manipulated using the Sirus 3-D viewing software set forth in the above-identified patents. 

[0043] To verify the accuracy of the 3-D imaging system, a string test phantom of known dimensions was imaged. 
The phantom is described in greater detail with reference to U.S. Patent No. 5,341,808, and was built from an acrylic 
frame with four horizontal planes 10mm apart. Surgical silk was stretched across each plane at spacings of 10mm. to 
produce a 10mm 3-D grid system. The phantom was submerged in a 7% glycerol solution (which mimics the speed of 

so sound in tissue) and scanned. Manual measurements obtained from the 3-D reconstructed image were compared with 
the known distances, and the results were found to be within ±0.1mm (or 1 %) for both the in-plane and between-plane 
readings. 

[0044] The algorithm was tested on wall-less carotid phantoms with varying degrees of stenoses in the internal 
carotid artery. These phantoms consisted of a "common" channel bifurcating into two additional channels, the "internal" 
55 and "external" carotid arteries, all of which were wall-less within a block of agar. The phantoms had stenoses of 30%. 
60%, and 70% as defined by NASCET (see above), and their channels were filled with a 7% glycerol solution prior to 
scanning. The phantoms were submerged in the glycerol solution so that about 2mm of the solution covered the sur- 
face, providing good acoustic coupling to the transducer. The axis of each phantom was manually adjusted to be 
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approximately parallel to the scan direction. Each phantom was scanned using B-mode imaging at a depth setting of 
2.8 cm and a single focal zone at 2.0 cm. 

[0045] After 3-D reconstruction of each set of acquired data, the vessel lumen was segmented. As discussed 
above, in the segmentation algorithm of the present invention, there are three user-specified parameters: the two force 

5 weightings w^, and w in , and the desired model resolution r. The internal and external weightings were set at w^, = 0.3 
and w irl = 1.0 in order to place greater emphasis on the internal contour smoothing, assuming the initialization process 
had performed well at localizing the edge boundaries. If greater weighting were to be given to the external gradient fea- 
tures, the contour could be attracted to high gradient noise in the image. The contour resolution was set at r = 15, giving 
an allowable range of distance between vertices of between 10 and 20 pixels. For the images analyzed, this corre- 

io sponded to a range of 0.8 to 1 7mm. 

[0046] A single seed paint in the center of any slice of the vessel is required to initiate the segmentation. To inves- 
tigate the variability in the resulting contour due to the placement of this seed, a slice was selected from the common 
carotid region and seeded using 15 different positions, ranging from the center to near the vessel boundary, as shown 
in Figure 6a. For each selected seed point, the initial contour was found and analyzed by determining the radial length 

ts from the approximate center of the lumen to the initial contour at 5 intervals. The results are plotted in Figure 6b. Fol- 
lowing each seed point selection, the segmentation proceeded as set forth above, resulting in the final contour. These 
contours were then analyzed as above for the initial contours, and the results are plotted in Figure 6c. In addition, the 
areas contained within the initial and final contours were determined, and their mean and standard deviation calculated. 
[0047] For each stenosed vessel, a single seed point was chosen at the approximate center of each branch. The 

20 segmentation procedure produced a contour every 0.25mm. The effects of the GDM application on a contour's energy 
is shown in Figure 7, while samples of the final segmentation are shown in Figures 8 and 9. The area within each con- 
tour was then calculated by summing the number of pixels within the contour and multiplying by the area of a pixel. 
These values were compared with the known areas from the numerical model used in casting the phantoms, and the 
results are shown in Figure 10 and Table II. 

25 



Table II 



Average magnitude difference between numerical and cal- 
culated areas. The stenosis severity is as defined by NAS- 
CET. 


Stenosis severity 


30% 


60% 


70% 


Common 


1.7% 


2.2% 


2.0% 


Internal 


1.6% 


2.2% 


2.0% 


External 


1.4% 


4.5% 


4.6% 


Total Average 


1.6% 


3.0% 


2.8% 



40 [0048] Because the curvature at a vertex in a GDM model is defined as the vector difference between the two 
edges joining a vertex, little or no curvature will be detected if the vertices are too closely spaced. Therefore, GDMs are 
limited in the resolution that they can provide. Accordingly, the effects on the final contour of varying the contour's res- 
olution was investigated using vertex spacings of r =8, 15, and 25 pixels. After convergence to a final contour, a plane 
above the bifurcation was selected and the results are shown in Figure 1 1. The areas within each contour were also 

45 calculated. 

[0049] In addition to measurement of accuracy, the time required to determine the initial contours for each vessel 
and the added time required to determine the final contour after application of the GDM, were measured. The timing 
was performed on a 9500 Power Macintosh prototype running at 100 MHz. 

[0050] Finally, in order to determine the algorithm's ability to segment images of human carotids, images from two 
so patients were analyzed. The images were acquired using two different acquisition techniques. One was scanned using 
the linear system described above m the phantom experiments, and the other was scanned using a freehand magnetic 
position and orientation measuring (POM) 3-D system (Rock of Birds, Ascension Technologies). The results of the seg- 
mentation are shown in Figure 12. 

[0051] The variability in segmentation caused by the seed point selection was evaluated using plots of the radial 
55 distance from the center of the structure to the contour, at five spacings. Figure 6b shows the radial contours resulting 
from the initialization steps. It can be seen that there are only small variations in the fifteen contours generated by the 
different seed points, even with respect to the extreme placements near the contour boundary. Figure 6c shows that the 
radial contours are smoother and more clustered after application of the GDM. The mean area of the contours was 



EP 1 030 187 A2 



47.7mm 2 after initialization and 47.1mm 2 after GDM application, showing that the area did not change significantly. 
However, the standard deviation of the 15 area estimates was reduced from 0.18mm 2 to 0.09mm 2 , showing that the 
application of the GDM resulted in the contours converging to a nearly common final result. 

[0052J Figure 7c shows a plot of the internal, external, and total force versus iteration number for a slice through the 

5 60% stenosed phantom at a location above the bifurcation. It can be seen that only five iterations are required for the 
deformation process to settle, likely due to the closeness of the initial contour. Also, it should be noted that the external 
force has increased slightly in the process of reducing the internal curvature of the contour. The corresponding initial 
and final contours are shown in Figures 7a and 7b, respectively. From these figures, it will be noted that the initial con- 
tour is a very close fit to the speckle at the lumen boundary, and that the GDM deformation serves mostly to smooth the 

io contour. Figure 8 shows the final segmentation for the 60% stenosed phantom in regions at the common carotid (8a), 
the bifurcation (8b), and above the bifurcation at the internal and external carotids (8c), while Figure 9 shows a cross- 
section through the length of the phantom. These figures provide a visual, qualitative assessment of the segmentation 
algorithm of the present invention, showing that the segmentation performs very well on the stenosed phantoms. 
[0053] Because the numerical model values for the cross-sectional areas were determined perpendicular to the 

is vessel axis (of the "common carotid"), an exact comparison with the scanned data would require alignment of the phan- 
tom with its axis parallel to the scanning axis. Although this was attempted visually, a small angular uncertainty 
remained. To minimize the effects of the angular uncertainty, separate comparisons were performed for each brunch of 
the three vessels. To match the scanned location with the numerical model, the measured cross-sectional areas were 
shifted and the chi-squared difference minimized. The results are plotted in Figure 10. From these figures, it is apparent 

20 that the cross-sectional areas match well with the numerical model. Quantitative comparisons between the measured 
and numerical areas are shown in Table II (above) for the brunches of the three vessels. From this table, it can be seen 
that the mean difference between the numerical and measured areas was 3% or less. 

[0054] The effects of the contour resolution in the GDM model on the final segmentation are shown in Figure 1 1 . 
As the between-vertex resolution increases (i.e. poorer resolution), more smoothing results as seen in the sequence 

25 from (a) to (d). However, if the resolution is allowed to be too large, the contour begins to lose its ability to closely follow 
the object boundary, as in (d). This is exemplified by the areas enclosed within the contours. For the left contour, the 
area decreased only slightly from 1 7.84 mm 2 for no GDM application (Figure 1 1 a) to 1 7.79 mm 2 for a vertex spacing of 
8 pixels (Figure 11 b), and to 17.34 mm 2 when a vertex spacing of 15 pixels (Figure 11 c) was applied. However, for a 
vertex spacing of 25 pixels (Figure 11d), an even larger decrease to 15.55 mm 2 occurred, as the spacings were too 

30 large and the model became too stiff for the contour to properly follow the boundary. 

[0055] The results for the segmentation time revealed that the determination of the initial contour for the entire ves- 
sel (200 slices) required only 21 seconds for each of the phantoms. To obtain a smoothed final contour for all slices in 
the phantom through the application of the GDM required only 4 additional seconds, due to the very few iterations 
(about 8) needed to bring the initial contours to a minimum energy configuration (Figure 7c). 

35 [0056] The avenge and standard deviation of the difference between the initial and final contours for each of the 
branches in each of the vessels is shown in Table III. 



Mean (A) and standard deviation (o(A)) of the difference in area (mm 2 ) between initial and final contours of each 
branch of each of the stenosed vessels. Branch 1 is the common, 2 the internal and 3 the external carotid. 


Stenosis severity 


30% 


60% 


70% 


Branch 


1 


2 


3 


1 


2 


3 


1 


2 


3 


A (mm') 


0.23 


0.39 


0.36 


0.23 


0.41 


0.34 


0.33 


0.20 


0.38 


a (A) (mm 2 ) 


0.31 


0.20 


0.30 


0.58 


0.23 


0.27 


0.26 


0.27 


0.21 



so [0057] The avenge difference was small, varying between 0.2 to 0.41 mm for each slice. This demonstrates that 
although the GDM smoothed the contour, it reduced the area by only a very small amount that was approximately equal 
to the standard deviation of the measurement. 

[0058] The results of the segmentation algorithm of the present invention applied to the images of the patient carot- 
ids are shown in Figures 12a and 1 2b, for the mechanically scanned and freehand acquired images, respectively. From 
55 these it can be concluded that the algorithm is capable of withstanding small areas of shadowing. Because the 2-D 
GDM used does not provide smoothing between contours, any fluctuations between acquired 2-D images will not be 
smoothed. This can be observed in both Figures 12a and 12b, where fluctuations due to cardiac pulsations in the com- 
mon carotid are present. 
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[0059] In summary, a fast automated segmentation method is provided according to the present invention with min- 
imal user interaction. The large variabilities that arise from user-defined initial contours have been eliminated, and 
experimental results have shown that the generated contours are insensitive to the placement within the vessel of the 
initial seed point needed to initiate the algorithm. The results of tests on carotid phantoms with a range of idealized sten- 

5 oses showed a difference of 3% or less from the numerical model used in making the phantoms. Exceptional speed was 
achieved in the initialization of the contour through the use of simple search techniques to identify a small number of 
candidate boundary points, which are further analyzed and compared using easily calculated features for each candi- 
date. Testing has shown that a total of only 25 seconds was needed for the analysis of a 200 slice, 23 MB volume on a 
9500 Power Macintosh prototype running at 100 MHz. 

10 [0060] Other embodiments and variations of the invention are contemplated. According to the current best mode, 
the candidates are compared through a simplistic weighting of their features. This technique works well for the images 
of phantoms and human carotids discussed herein above. However, in order to extend the use of the algorithm to a 
wider range of images of differing quality, additional edge features may be required, as well as a more explicit knowl- 
edge-based analysis of feature combinations. A neural network is well suited for this application, as it allows for a non- 

15 linear weighting of the features to be learned. Additionally, the neural network could simultaneously be presented with 
data from previously selected candidates in adjacent rays, both in the current slice as well as adjacent slices, giving it 
a larger and more regional view in the learning and scoring of boundary candidates. The initial contour generation pres- 
ently uses data from tree slices, while the GDM minimizes the contour's energy with data in the current slice only, with- 
out regard to continuity in the third dimension. To maintain and enhance the limited 3-D continuity of the initial contours, 

20 the GDM model can be extended to 3-D. This would provide increased robustness, enabling the contour to better over- 
come large areas of shadow, speckle, and reflection, all of which are problems frequently encountered in imaging 
human carotids. Finally, with successful segmentation and the availability of the cross-sectional areas along the length 
of the vessel, the next step is to use this information to arrive at a quantitative measurement of stenosis. Currently, x- 
ray angiography is the method of choice for quantitative assessment. With angiography, stenosis is described as the 

25 ratio of the diameter at the point of greatest stenosis to that at a location distal to the stenosis. This measure has been 
found to be a good indicator for identifying those at increased risk for stroke. In x-ray imaging, the magnification factors 
of the imaged structure are widely varying and unknown, making it impossible to directly measure lumen size. One-rea- 
son for taking the ratio Of the two diameters is to cancel out the effects of the magnification. With ultrasound, it is difficult 
to obtain the same ratio measurement, as the jaw bone limits the viewing length along the carotid arteries. However, 

30 information on voxel dimensions is available and can be used in calculating the true lumen areas. The measurement for 
stenosis can thus simply be taken as the cross-sectional area at the point of greatest stenosis. The availability of such 
quantitative information can improve the diagnostic capabilities of ultrasound imaging.allowing it the potential for play- 
ing an increasing role in the monitoring of carotid disease progression. 

[0061] A person skilled in the art of 3-D medical ultrasound imaging may conceive of other variations and modifica- 
35 Sons of the invention, all of which are believed to be within the sphere and scope of the invention as set forth in the 
appended claims. 

Claims 

40 1. An automated segmentation method for generating a surface contour from an ultrasound image of a fluid filled 
region, comprising the steps of: 

a) identifying an initial contour by 

45 specifying a single seed point approximately centrally in said ultrasound image of said fluid filled region. 

extending a plurality of radial lines from said single seed point, each of said radial lines being characterized 
by an intensity profile derived from said image, 

performing a threshold analysis of said intensity profile for each of said radial lines and in response select- 
ing a plurality of initial boundary point candidates along respective ones of said radial lines, 
so comparing respective ones of said plurality of initial boundary point candidates along respective ones of 

said radial lines with further initial boundary point candidates along adjacent radial lines and selecting from 
among said plurality of initial boundary point candidates one boundary point per radial line which is con- 
tinuous with said further initial boundary point candidates along said adjacent radial lines, whereby said 
initial contour is identified as a set of selected boundary points connected by respective edges; and 

55 

b) applying a geometrically deformable model to said initial contour for smoothing said edges and thereby gen- 
erating said surface contour. 
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2. The method of claim 1, wherein said step of comparing respective ones of said plurality of initial boundary point 
candidates along respective ones of said radial lines is conducted with respect to two adjacent radial lines in a sin- 
gle 2D plane of said image and with respect to three adjacent radial lines in each adjacent plane, such that a con- 
tinuity constraint is imposed on selection of said one boundary point per radial line in three dimensions. 

5 

3. The method of claim 2, wherein individual ones of said plurality of radial lines extending from said single seed point 
are separated by 5°. 

4. The method of claim 1 , wherein said step of performing a threshold analysis of said intensity profile for each of said 
10 radial lines further includes evaluating said intensity profile at a plurality of thresholds and in response generating 

a multilevel binary wee of values wherein a "0" indicates a voxel intensity below one of said thresholds and a "1" 
indicates a voxel intensity greater than one of said thresholds. 

5. The method of claim 4, wherein said thresholds are automatically selected for each one of said radial lines, thereby 
is rendering said thresholds adaptive to variable properties of said image. 

6. The method of claim 5, wherein said thresholds are calculated by identifying a voxel on said radial line character- 
ized by maximum intensity, l m , as occupying at least a portion of the initial boundary along said radial line, and 
selecting said thresholds such that said voxel occupies at least four levels of said multilevel binary tree. 

20 

7. The method of claim 6. wherein said thresholds T, are calculated using an intensity value l 0 at said seed point and 
said maximum intensity l m as follows: 

25 

T, = l 0 + (i*AI);i = 1,2,3,4.5. 

8. The method of claim 7, further comprising a step of increasing Al in the event of an excessive number of said 
boundary point candidates, according to Al = (I m - 1 0 )/3 , and thereafter reevaluating said binary tree. 

30 

9. The method of claim 8. further comprising the steps of determining start position, length and number of levels of 
said binary tree occupied by said voxel, rating said voxel according to four criteria which are (i) edge position rela- 
tive to said further additional boundary point candidates, (ii) number of levels of said binary tree occupied, (iii) edge 
length, and (iv) variance in the voxels along said radial line between the seed point and boundary point candidate, 

35 and in response selecting said one boundary point. 

10. The method of claim 9, wherein, within a set of candidates in a radial line, each criterion is normalized across all 
candidates to give a value between 0 and 1, for ensuring that at least one candidate receives a maximum value in 
each of the criteria. 

AO 

1 1 . An apparatus for generating a surface contour from an ultrasound image of a fluid filled region by automated seg- 
mentation, comprising: 

means for obtaining an ultrasound image; 
45 processing means for identifying an initial contour based on a specified single seed point approximately cen- 

trally in said ultrasound image of said fluid filled region, by 

extending a plurality of radial lines from said single seed point, each of said radial lines being characterised by 
an intensity profile derived from said image, 

performing a threshold analysis of said intensity profile for each of said radial lines and in response selecting 
so a plurality of initial boundary point candidates along respective ones of said radial lines, 

comparing respective ones of said plurality of initial boundary point candidates along respective ones of said 
radial lines with further initial boundary point candidates along adjacent radial lines and selecting from among 
said plurality of initial boundary point candidates one boundary point per radial line which is continuous with 
said further initial boundary point candidates along said adjacent radial lines, whereby said initial contour is 
55 identified as a set of selected boundary points connected by respective edges; 

processing means for applying a geometrically deformable model to said initial contour for smoothing said 
edges and thereby generating said surface contour; and 
means for displaying said generated surface contour. 
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